clear all  %#ok

Consolidate = 0 ;
if Consolidate==1
    diary(strcat('.\BOOTSTRAPS\CONSILIDATEandPROCESS_BOOTSTRAPS_LOGMSM7.txt'))
    load('PointEstimates_MSM7.mat')
    BootstrapResults = cell(1,1) ;
    format long
    LoadVarious = 1 ;
    if LoadVarious>0 && LoadVarious<=1.5
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%  BELOW IS EXAMPLE CODE OF HOW DIFFERENT BOOTSTRAP BATCH FILES CAN BE CONSOLIDATED.  WE DID
%%%%%%%  NOT INCLUDE THE BATCH FILES BECAUSE OF THEIR LARGE SIZE (29GB TOTAL).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        tic
        disp('loading the first 100 bootstrap results...')
        load('.\BOOTSTRAPS\BSResultsMSM7_1_100.mat');
        for s=1:100; BootstrapResults{s,1} = BSResultsMSM_1_100{s,1}; end
            clear BSResultsMSM_1_100
        toc
        tic
        disp('loading 101...200 bootstrap results...')
        load('.\BOOTSTRAPS\BSResultsMSM7_101_200.mat');
        for s=101:200; BootstrapResults{s,1} = BSResultsMSM_101_200{s,1}; end
            clear BSResultsMSM_101_200
        toc
        tic
        disp('loading 201...300 bootstrap results...')
        load('.\BOOTSTRAPS\BSResultsMSM7_201_300.mat');
        for s=201:300; BootstrapResults{s,1} = BSResultsMSM_201_300{s,1}; end
            clear BSResultsMSM_201_300
        toc
        tic
        disp('loading 301...400 bootstrap results...')
        load('.\BOOTSTRAPS\BSResultsMSM7_301_400.mat');
        for s=301:400; BootstrapResults{s,1} = BSResultsMSM_301_400{s,1}; end
            clear BSResultsMSM_301_400
        toc
        tic
        disp('loading 401...450 bootstrap results...')
        load('.\BOOTSTRAPS\BSResultsMSM7_401_450.mat');
        for s=401:450; BootstrapResults{s,1} = BSResultsMSM_401_450{s,1}; end
            clear BSResultsMSM_401_450
        toc
        
        

        toc/60 %#ok
        clear s BSResults*
        BSResults = BootstrapResults ;
        clear BootstrapResults
        for s=1:length(BSResults)
            BSResults{s,1}.BSSampID = s ;
        end
        %%%%%%  ALTER THE FOLLOWING LINE TO KEEP THE FIRST 400 OF THE 450 BOOTSTRAPS ABOVE THAT
        %%%%%%  TERMINATED.  OCCASIONALLY A BOOTSTRAP SAMPLE WILL NOT TERMINATE OR WILL PRODUCE
        %%%%%%  NAN'S IN THE OUTPUT.  THIS IS WHY WE RUN 450 TO GET A FINAL SAMPLE OF 400
        %%%%%%  BOOTSTRAPS. 
        if LoadVarious==1
            msg = '**************************PROCESSBOOTSTRAPS INSTRUCTIONS: Step 1 consolidation complete.  Now, check BSResults in the workspace to determine what set of bootstrap indices gives the first 450 bootstrap estimates across all batches run, and update the right-hand side of the next line of code below accordingly.  When you are done, re-set the value of LoadVarious to 1.5 and re-run ProcessBootstraps.' ;
            error('ProcessBootstraps:Step1',msg)
        elseif LoadVarious==1.5
            BSResults = BSResults(1:400,1);
            save('.\BOOTSTRAPS\BSResultsMSM7_1_450.mat','BSResults','-v7.3')
            msg = '**************************PROCESSBOOTSTRAPS INSTRUCTIONS: Step 2 consolidation complete.  Now, re-set the value of LoadVarious to 2 and re-run ProcessBootstraps.' ;
            error('ProcessBootstraps:Step2',msg)
        end
    elseif LoadVarious==2
        load('.\BOOTSTRAPS\BSResultsMSM7_1_450.mat');
        disp(strcat('loading_',num2str(length(BSResults)),'_bootstrap results...'))
        for s=1:length(BSResults)
            BSResults{s,1}.BSSampID = s ; %#ok
        end

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%% THE FOLLOWING BLOCK OF CODE PERFORMS ONE LAST SANITY CHECK: IT CHECKS THE SOLUTION
        %%%%%%% OF EVERY BOOTSTRAP AGAINST ALL OTHER BOOTSTRAP SOLUTIONS TO SEE IF AN IMPROVEMENT
        %%%%%%% CAN BE FOUND.  IF SO, THEN IT RUNS ONE LAST ROUND OF PATTERN SEARCH FROM THE
        %%%%%%% IMPROVED STARTING POINT.
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        SanityCheck = 1 ;
        diary(strcat('.\BOOTSTRAPS\BSResultsMSM7_1_400_Sanity_LOG.txt'))
        load('PointEstimates_MSM7.mat')
        if SanityCheck==1
            Solutions = cell(length(BSResults),1) ;
            for tt=1:length(BSResults)  %
                temp = size(fieldnames(BSResults{tt,1})) ;
                if temp(1)==3
                    Solutions{tt,1} = BSResults{tt,1}.LPointEst.params ;
                end
            end


            start  = 1 ;
            finish = length(Solutions(:,1)) ;
            BestSol   = cell(length(Solutions(:,1)),1) ;
            for tt = 1:length(Solutions(:,1))
                if ~isempty(Solutions{tt,1})
                    tic
                    disp(strcat('Now checking solution_',num2str(tt),'_of_',num2str(length(Solutions(:,1))),'_bootstrap estimates (#',num2str(tt),'_of_',num2str(1095),')...'))
                    datetime(now,'ConvertFrom','datenum') %#ok

                    ssrlog = nan(length(Solutions(:,1)),1) ;
                    student_id   = BSResults{tt,1}.LPointEst.student_id ;
                    Q            = BSResults{tt,1}.LPointEst.Q ;
                    T            = BSResults{tt,1}.LPointEst.T ;
                    ThetaE       = BSResults{tt,1}.LPointEst.ThetaE ;
                    contract     = BSResults{tt,1}.LPointEst.contract ;
                    Q80freq      = BSResults{tt,1}.LPointEst.Q80freq ;
                    pay          = BSResults{tt,1}.LPointEst.pay ;
                    SimTimeCum   = BSResults{tt,1}.LPointEst.SimTimeCum ;
                    Activeids    = BSResults{tt,1}.LPointEst.Activeids ;
                    Marginalids1 = BSResults{tt,1}.LPointEst.Marginalids1 ;
                    Marginalids2 = BSResults{tt,1}.LPointEst.Marginalids2 ;
                    Tdom         = BSResults{tt,1}.LPointEst.Tdom ;
                    EQgivenT_lb  = BSResults{tt,1}.LPointEst.EQgivenT_lb ;
                    EQgivenT_ub  = BSResults{tt,1}.LPointEst.EQgivenT_ub ;
                    cutoffs      = BSResults{tt,1}.EPointEst.cutoffs ;
                    BSstate      = [1,max(LPointEst.QuantileFunq1),max(LPointEst.QuantileFunq2),max(LPointEst.QuantileFunq3)] ;
                    tic
                    parfor ss=1:length(Solutions(:,1))
                        try
                            ssrtemp = ObjFun_cSGMM(Solutions{ss,1},LPointEst.knotvec_c,student_id,Q,T,...  %%%
                                ThetaE,contract,Q80freq,pay,SimTimeCum,... %%%
                                Activeids,Marginalids1,Marginalids2,...
                                Tdom,EQgivenT_lb,EQgivenT_ub,cutoffs,... %%%
                                [],[],[],[],0,BSstate) ;
                            ssrlog(ss)  = ssrtemp ;
                        catch
                            % % % % %                     disp('!!!')
                            ssrlog(ss) = inf ;
                        end
                    end
                    toc
                    if ssrlog(tt) == min(ssrlog(start:finish))
                        BestSol{tt} = tt ;
                        disp(strcat('Bootstrap_',num2str(tt-start+1),'_of_',num2str(finish-start+1),'_is *********OPTIMAL---OPTIMAL---OPTIMAL********* (#',num2str(tt),'_of_',num2str(BootstrapSims),')...'))

                    else
                        BestSol{tt} = find(ssrlog==min(ssrlog(start:finish))) ;
                        disp(strcat('Bootstrap_',num2str(tt-start+1),'_of_',num2str(finish-start+1),'_is !!!!!!!!!!!!!!!!!!SUB-OPTIMAL!!!!!!!!!!!!!!!!!! (#',num2str(tt),'_of_',num2str(BootstrapSims),')...'))
                        format long
                        pi_0 = Solutions{tt,1} ;
                        pi_better = Solutions{BestSol{tt},1} ;
                        CheckAltGuesses = [pi_0 pi_better 100*(pi_better-pi_0)./pi_0;...
                            ssrlog(tt) ssrlog(BestSol{tt}) 100*(ssrlog(BestSol{tt})-ssrlog(tt))/ssrlog(tt)]   %#ok


                        BSResults{tt} = struct ; %#ok
                        sampindx     = [] ;
                        for kk=1:length(BSSamp{tt,1})
                            temp     = find(ProductionDataNonCum.sid==BSSamp{tt,1}(kk)) ;  %%%%This is all the rows corresponding to the current re-sampled test subject
                            sampindx = [sampindx; [temp, kk*ones(size(temp))]] ; %#ok  %%%%The second column will be a new student id index to replace the old one.  This is necessary since subjects are re-sampled with replacement and therefore can be selected more than once.  Thus, each re-sampled copy of the same student requires a unique new index.
                        end
                        ProductionSamp     = ProductionDataNonCum(sampindx(:,1),:) ;
                        ProductionSamp.sid = sampindx(:,2) ;  %%%%This replaces the old indices with the new indexing for the re-sampled collection of students.
                        
                        EPointEsttemp                     = ThetaEEstimator(ProductionSamp,0) ;
                        BSResults{tt}.EPointEst           = EPointEsttemp ; %#ok
                        cutoffs                           = EPointEsttemp.cutoffs ; %#ok

                        sampindx = nan(length(BSSamp{tt,1}),2) ;
                        for kk=1:length(BSSamp{tt,1})
                            sampindx(kk,1) = find(AchievementData.id==BSSamp{tt,1}(kk)) ;  %%%%This is all the rows corresponding to the current re-sampled test subject
                            sampindx(kk,2) = kk ;  %%%%The second column will be a new student id index to replace the old one.  This is necessary since subjects are re-sampled with replacement and therefore can be selected more than once.  Thus, each re-sampled copy of the same student requires a unique new index.
                        end
                        AchievementSamp           = AchievementData(sampindx(:,1),:) ;
                        AchievementSamp.id        = sampindx(:,2) ;  %%%%This replaces the old indices with the new indexing for the re-sampled collection of students.
                        AchievementSamp.ThetaE(:) = NaN ;
                        for kk=1:length(BSResults{tt}.EPointEst.studlistFULL)
                            AchievementSamp.ThetaE(AchievementSamp.id==BSResults{tt}.EPointEst.studlistFULL(kk)) = BSResults{tt}.EPointEst.ThetaE(kk) ;
                        end
                        
                        knotvec_c = LPointEst.knotvec_c; %#ok
                        disp('%%%%PSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPS%%%%%')
                        disp(strcat('%%%% TERMINAL ROUND OF SIMULATED GMM ESTIMATION VIA PATTERN SEARCH...     %%%%%'))
                        disp('%%%%PSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPS%%%%%')
                        BSstate(1) = 3 ;
                        LPointEsttemp = ThetaLEstimatorSGMM(EPointEsttemp,AchievementSamp,LPointEst.knotvec_c,[],pi_better,0,BSstate,0,[]) ;
                        ssrPSterminal = LPointEsttemp.ssrvalc ; %#ok
                        format long
                        LPointEsttemp.params
                        BSResults{tt}.LPointEst = LPointEsttemp ; %#ok
                    end
                    toc
                end
            end
        end
        flag = zeros(length(BSResults),1); 
        for ii=1:length(BSResults) 
            if isempty(Solutions{ii,1}) 
                flag(ii)=1; 
            end
        end
        keepindx = find(flag==0) ;
        keepindx(keepindx>400) = [] ;  %%%%%We only need to keep the first 400 bootstraps.  Having a number of bootstraps that is divisible by 4 simplifies computation of confidence intervals later...
        BSResults = BSResults(keepindx) ;  %%%%Here we drop the small number of instances where the solver couldn't find a solution
        save('.\BOOTSTRAPS\BSResultsMSM7_1_400.mat','BSResults','-v7.3')
        msg = '**************************PROCESSBOOTSTRAPS INSTRUCTIONS: Step 3 consolidation complete.  Now, re-set the value of LoadVarious to 0 and re-run ProcessBootstraps.  This will execute a The next phase where we compute productivity and motivation fixed effects for active students.' ;
        error('ProcessBootstraps:Step3',msg)
    else 
        load('.\BOOTSTRAPS\BSResultsMSM7_1_400.mat')
    end


    
    %%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%THIS SEGMENT OF CODE PRODUCES ESTIMATES OF STUDENT FIXED EFFECTS (PRODUCTIVITY AND 
    %%%%%%MOTIVATION) FOR EACH OF THE BOOTSTRAP COST FUNCTIONS COMPUTED AND PROCESSED ABOVE.
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    diary(strcat('.\BOOTSTRAPS\TypeBootstrap_LOGMSM7.txt'))
    knotvec_c   = LPointEst.knotvec_c;
    pi_c        = LPointEst.params ;
    NBS         = length(BSResults) ;
    signif      = 0.05 ; %#ok
    BSestimates = nan(length(pi_c)+3,NBS) ;

    for s=1:NBS
        BSestimates(1:length(pi_c),s)   = BSResults{s,1}.LPointEst.params ;
        BSestimates(length(pi_c)+1,s)   = BSResults{s,1}.EPointEst.tau0 ;
        BSestimates(length(pi_c)+2,s)   = BSResults{s,1}.EPointEst.tau1 ;
        BSestimates(length(pi_c)+3,s)   = BSResults{s,1}.EPointEst.phi ;
    end


    BSestimates(:,isnan(BSestimates(1,:))) = [] ; %#ok


    

    %%%%Point Values of ThetaL for full-output students (Q==80) will come from the Tobit Results.  For those students, this
    %%%%variable will contain an upper bound on ThetaL (i.e., assuming the actual stop point would have been exactly Q==80).
    ThetaLBS         = nan(length(LPointEst.ThetaE),NBS) ;
    ThetaEBS         = nan(length(LPointEst.ThetaE),NBS) ;
    Nstud            = length(LPointEst.student_id) ;
    Q                = LPointEst.Q ;
    T                = LPointEst.T ; %#ok
    student_id       = LPointEst.student_id ;
    contract         = LPointEst.contract ;
    ThetaEThetaLBSid =  student_id ; %#ok

    id  = ProductionDataNonCum.sid ;
    tau = ProductionDataNonCum.t ;
    q   = ProductionDataNonCum.k ;
    dropindx = [] ;
    for s=1:NBS
        disp('***')
        disp('***')
        disp('***')
        disp(strcat('Now doing types post-estimation for Bootstrap #',num2str(s),'...'))
        pi_cBS       = BSResults{s,1}.LPointEst.params' ;
        tau0temp         = BSResults{s,1}.EPointEst.tau0 ;
        tau1temp         = BSResults{s,1}.EPointEst.tau1 ;
        phitemp          = BSResults{s,1}.EPointEst.phi ;

        if EPointEst.BurninModel==1
            burnin = 1*(q==1) ;
            Z = log(tau) - tau1temp*burnin + phitemp*log(q) ;
            for ii=1:length(student_id)
                if Q(ii)>=2
                    ThetaEBS(ii,s) = exp( sum( Z(id==student_id(ii)) - log(tau0temp) )/Q(ii) ) ;
                end
            end
        end

        disp('computing preliminaries...')
        tic
        P = ThetaLEstimatorSGMM(BSResults{s,1}.EPointEst,AchievementData,...
            knotvec_c,[],pi_cBS,0,[4;max(LPointEst.Q80freq(:,1,1));max(LPointEst.Q80freq(:,1,2));max(LPointEst.Q80freq(:,1,3))],1) ;
        toc
        disp('now computing Theta_l types...')
        tic
        try
            [~,out] = ObjFun_cSGMM(pi_cBS,knotvec_c,P.student_id,P.Q,P.T,...  %%%
                P.ThetaE,P.contract,P.Q80freq,P.pay,P.SimTimeCum,... %%%
                P.Activeids,P.Marginalids1,P.Marginalids2,...
                P.Tdom,P.EQgivenT_lb,P.EQgivenT_ub,P.cutoffs,... %%%
                [],[],[],[],0,[4;max(LPointEst.Q80freq(:,1,1));max(LPointEst.Q80freq(:,1,2));max(LPointEst.Q80freq(:,1,3))]) ;
            toc
            for ii=1:Nstud
                if Q(ii)<80 && Q(ii)>1
                    ThetaLBS(ii,s) = out.ThetaLACT(out.student_id==student_id(ii)) ;
                elseif Q(ii)==1
                    if contract(ii)==1
                        ThetaLBS(ii,s) = out.ThetaLMRG1( find(out.student_id==student_id(ii)) - length(out.ThetaLACT) ) ;
                    elseif contract(ii)==2
                        ThetaLBS(ii,s) = out.ThetaLMRG1( find(out.student_id==student_id(ii)) - length(out.ThetaLACT) - length(out.ThetaLMRG1) ) ;
                    end
                elseif Q(ii)==80
                    ThetaLBS(ii,s) = out.ThetaL80(out.ThetaL80(:,1)==student_id(ii),2) ;
                end
            end
        catch
            toc
            disp(strcat('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  Error occurred in bootstrap number_',num2str(s),'.  Skipping that one for the purpose of ThetaL computation...'))
            dropindx = [dropindx; s] ; %#ok
        end
        
    end
    clear ii id tau q paramstemp Q80freq UdomHettemp F_usmoothHettemp
    clear cutoffstemp tau1temp phitemp Z indx ThetaL80temp1 ThetaL80temp2 ThetaL80temp3
    clear ThetaL80temp4 ThetaL80temp5 Ltemp L80temp s Nstud

    format shortg

    save('.\BOOTSTRAPS\ProcessBootstrapsMSM7.mat','-v7.3')
    %%
    diary off
    msg = '**************************PROCESSBOOTSTRAPS INSTRUCTIONS: Step 4 consolidation complete: BOOTSTRAPS ARE NOW CONSOLIDATED, AND PRODUCTIVITY MOTIVATION TYPES ARE COMPUTED FOR EACH ACTIVE STUDENT IN THE ORIGINAL SAMPLE, GIVEN EACH BOOTSTRAPPED ESTIMATE OF COMMON PARAMETERS.  Now, re-set the value of LoadVarious to 1, re-set the value of Consolidate to 0, and re-run ProcessBootstraps.  This will move to the next phase where confidence intervals are computed for common parameters, and some plots are produced.' ;
    error('ProcessBootstraps:Step4',msg)
else
    load('.\BOOTSTRAPS\ProcessBootstrapsMSM7.mat')
    load('g2.mat')
    clear BSestimatesBC
%     vvvvvvvvv
    dropindx(dropindx(:,1)>length(BSResults)) = [] ;
    extraneeded = sum(dropindx<=400) ;
    if sum(dropindx<=400+extraneeded) == extraneeded
        dropindx = unique([dropindx; (400+extraneeded+1:length(BSResults))']) ;
    end
    keepset = (1:length(BSResults))' ;
    keepset(dropindx) = [] ;
    dropend = mod(length(keepset),40) ; %%%Keeping total sample size at multiples of 40 simplifies computation of the 95\% bounds
    keepset = keepset(1:end-dropend) ;
    BSestimates(:,dropindx) = [] ; %#ok
    BSResults               = BSResults(keepset,1) ;
    NBS                     = length(keepset) ;
    ThetaEBS(:,dropindx)    = [] ;
    ThetaLBS(:,dropindx)    = [] ;
    BSSampkeepset           = nan(NBS,1) ;
    for ii=1:NBS
        BSSampkeepset(ii,1) = BSResults{ii,1}.BSSampID ;
    end
    BSSamp        = BSSamp(BSSampkeepset,1) ;
    BootstrapSims = length(BSResults) ;
    %%%%The following loop is an attempt to clear memory space:
    for ii=1:NBS
        BSResults{ii,1}.EPointEst = rmfield(BSResults{ii,1}.EPointEst,{'U','ThetaEwithU','FirstTime','UdomHet','F_usmoothHet','ShockDistMassMarginal'}) ;
        BSResults{ii,1}.LPointEst = rmfield(BSResults{ii,1}.LPointEst,{'domq1','QuantileFunq1','domq2','QuantileFunq2','domq3','QuantileFunq3','Q80freq','outputcpat2','ssrvalc',...
                                                                       'ssr_1lowtail','ssr_1mid','ssr_1upptail','guardrail1up','guardrail1lo',...
                                                                       'ssr_2lowtail','ssr_2mid','ssr_2upptail','guardrail2up','guardrail2lo',...
                                                                       'ssr_3lowtail','ssr_3mid','ssr_3upptail','guardrail3up','guardrail3lo',...
                                                                       'ssr'}) ;
    end


    Nstud            = length(LPointEst.student_id) ;
    Q                = LPointEst.Q ;
    T                = LPointEst.T ;
    student_id       = LPointEst.student_id ;
    contract         = LPointEst.contract ;
    ThetaEThetaLBSid =  student_id ;
    medThetaL = median(LPointEst.ThetaLACT(:,1),'omitnan') ;
    NBS       = length(BSResults) ;


    disp('normalizing all stage-2 Bootstrap estimates to have the same median ThetaL type for active students...')
    ThetaLBS_norm = nan(size(ThetaLBS)) ;
    for ss=1:NBS
        BSResults{ss,1}.LPointEst.params_norm = BSResults{ss,1}.LPointEst.params*(median(ThetaLBS(Q>=2,ss))/medThetaL) ;
        ThetaLBS_norm(:,ss)                   = ThetaLBS(:,ss)*(medThetaL/median(ThetaLBS(Q>=2,ss))) ;
    end

    disp('computing confidence intervals for stage-1 and stage-2 parameters...')
    signif      = 0.05 ; %#ok
    plotsignif  = 0.1 ;
    BSestimates = nan(length(pi_c)+3,NBS) ;
    for s=1:NBS
        BSestimates(1:length(pi_c),s) = BSResults{s,1}.LPointEst.params ;
        BSestimates(length(pi_c)+1,s) = BSResults{s,1}.EPointEst.tau0 ;
        BSestimates(length(pi_c)+2,s) = BSResults{s,1}.EPointEst.tau1 ;
        BSestimates(length(pi_c)+3,s) = BSResults{s,1}.EPointEst.phi ;
    end
    BSestimates(:,isnan(BSestimates(1,:))) = [] ;
    BSestimatesS                           = sort(BSestimates,2) ; % - mean(BSestimates,2) + [LPointEst.params;EPointEst.tau0;EPointEst.tau1;EPointEst.phi] ;
    NBS             = length(BSestimates(1,:)) ;
    format shortg
    signif = 0.9 ;
    CImethod0       = [[nan(2,1); knotvec_c; nan(3,1)] BSestimatesS(:,round(NBS*(signif/2))) [LPointEst.params; EPointEst.tau0; EPointEst.tau1; EPointEst.phi] BSestimatesS(:,round(NBS*(1-signif/2)))] ;
    BSestimatesBC = sort(BSestimates - median(BSestimates,2) + [LPointEst.params; EPointEst.tau0; EPointEst.tau1; EPointEst.phi],2) ;  %%%%This is the bias-corrected sample of B-spline parameters
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%  
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    CImethod1       = [[nan(2,1); knotvec_c; nan(3,1)] BSestimatesBC(:,round(NBS*(signif/2))) [LPointEst.params; EPointEst.tau0; EPointEst.tau1; EPointEst.phi] BSestimatesBC(:,round(NBS*(1-signif/2)))]  %#ok
    if length(knotvec_c)==8  
        StructuralParams={'Knots','VARS', 'Lower95CI','PointEst','UPPer95CI';... 
                          CImethod1(1,1), 'pi_c1', CImethod1(1,2), CImethod1(1,3), CImethod1(1,4);...
                          CImethod1(2,1), 'pi_c2', CImethod1(2,2), CImethod1(2,3), CImethod1(2,4);...
                          CImethod1(3,1), 'pi_c3', CImethod1(3,2), CImethod1(3,3), CImethod1(3,4);...
                          CImethod1(4,1), 'pi_c4', CImethod1(4,2), CImethod1(4,3), CImethod1(4,4);...
                          CImethod1(5,1), 'pi_c5', CImethod1(5,2), CImethod1(5,3), CImethod1(5,4);...
                          CImethod1(6,1), 'pi_c6', CImethod1(6,2), CImethod1(6,3), CImethod1(6,4);...
                          CImethod1(7,1), 'pi_c7', CImethod1(7,2), CImethod1(7,3), CImethod1(7,4);...
                          CImethod1(8,1), 'pi_c8', CImethod1(8,2), CImethod1(8,3), CImethod1(8,4);...
                          CImethod1(9,1), 'pi_c9', CImethod1(9,2), CImethod1(9,3), CImethod1(9,4);...
                          CImethod1(10,1), 'pi_c10', CImethod1(10,2), CImethod1(10,3), CImethod1(10,4);...
                          CImethod1(11,1), 'tau0', CImethod1(11,2), CImethod1(11,3), CImethod1(11,4);...
                          CImethod1(12,1), 'tau1', CImethod1(12,2), CImethod1(12,3), CImethod1(12,4);...
                          CImethod1(13,1), 'phi', CImethod1(13,2), CImethod1(13,3), CImethod1(13,4)}  %#ok
                      
    end

    %%
    rowindx = find(student_id==1511) ;  %%%%kid with median ThetaL type
    testdom   = linspace(knotvec_c(1),knotvec_c(end),2000)' ;
    testdom   = unique([testdom; 45;90;180;360;720]) ;
    Cmat      = nan(length(testdom),NBS) ;

    for s=1:NBS
        Ctemp = BsplineEval3(knotvec_c,BSestimates(1:length(knotvec_c)+2,s),testdom) ;
        Cmat(:,s)      = ThetaLBS(rowindx,s)*Ctemp ;

    end
    C      = BsplineEval3(knotvec_c,pi_c,testdom) ;
    C               = LPointEst.ThetaLACT(rowindx)*C ;

    medianCmat      = median(Cmat,2,'omitnan') ;

    CmatSORT        = nan(size(Cmat)) ;

    for s=1:NBS
        CmatSORT(:,s)      = Cmat(:,s) - medianCmat + C ;

    end
    CmatSORT      = sort(CmatSORT,2) ;

    nbs         = sum(~isnan(CmatSORT(1,:))) ;
    scale = 1 ; %median(AchievementData.ThetaL(AchievementData.Q>=2)) ;





    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%                 FIGURE 2:                 %%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    figure;
    hold on
    plot(testdom,(scale*C),'b','LineWidth',3);
    plot(testdom,(scale*CmatSORT(:,round(nbs*(plotsignif/2))+1)),'--k')
    plot(testdom,(scale*CmatSORT(:,round(nbs*(1-plotsignif/2))-1)),'--k')
    scatter(T(rowindx),interp1(testdom,(scale*C),T(rowindx)),150,'o','filled','MarkerEdgeColor','k','MarkerFaceColor','r','LineWidth',2)
    xlim([0,700])
    ylim([0,150])
    xlabel('10-DAY TOTAL (EXTRACURRICULAR) WORK TIME t (minutes)','FontWeight','bold','FontSize',20)
    ylabel({'COST LEVELS, C(t;\theta_{mi})';'(dollars, scaled to median active \theta_{mi})'},'FontWeight','bold','FontSize',20)
    legend('POINT ESTIMATE',strcat(['UPPER',' ',num2str(round(100*(1-plotsignif))),'% CONF. BAND (for median active \theta_{mi})']),...
                            strcat(['LOWER',' ',num2str(round(100*(1-plotsignif))),'% CONF. BAND (for median active \theta_{mi})']),...
                            strcat(['=',num2str(round(T(rowindx),1)),',',' ','OBSERVED STOP TIME (for median active \theta_{mi})']),'Location','SouthEast','FontWeight','bold','FontSize',15)
    box on
    grid on

    xref = [117; 65];
    markerpoints = [150;200;300;450;600];
    for ii=1:5
        xval = (xref(1)+markerpoints(ii))/(xref(1)+720+xref(2)) ;
        x = [xval xval];
        if ii==1
            y = [0.335 0.19];
        elseif ii==2
            y = [0.77 0.33];
        elseif ii==3
            y = [0.455 0.305];
        elseif ii==4
            y = [0.605 0.455];
        elseif ii==5
            y = [0.82 0.67];
        end
        if ii~=2 && ii<5
            a = annotation('textarrow',x,y,'String',strcat('C(',num2str(markerpoints(ii)/10),' min/day)=$',num2str(round(scale*interp1(testdom,C,markerpoints(ii),'linear')/10,2)),'/day')) ;
            a.FontSize = 16 ;
            a.FontWeight = 'bold' ;
        elseif ii==5
            a = annotation('textarrow',x,y,'String',strcat('C(1 hr/day)=$',num2str(round(scale*interp1(testdom,C,markerpoints(ii),'linear')/10,2)),'/day')) ;
            a.FontSize = 16 ;
            a.FontWeight = 'bold' ;
        end
    end

    %%




    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%            FIGURE 3 (LEFT PANEL):         %%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    figure;
    hold on
    scale = 1 ;
    plot(testdom,log10(scale*C),'r','LineWidth',3);
    x = [0.47 0.3] ;
    y = [0.305 0.59] ;
    a = annotation('textarrow',x,y,'String',strcat('C(15 min/day)=$',num2str(round(scale*interp1(testdom,C,150,'linear')/10,2)),'/day')) ;
    a.FontSize = 20 ;
    a.FontWeight = 'bold' ;
    x = [0.67 0.7925] ;
    y = [0.5 0.75] ;
    a = annotation('textarrow',x,y,'String',strcat('C(1 hr/day)=$',num2str(round(scale*interp1(testdom,C,600,'linear')/10,2)),'/day')) ;
    a.FontSize = 20 ;
    a.FontWeight = 'bold' ;
    [Ftemp,Xtemp] = ecdf(AchievementData.ThetaL(AchievementData.Q>=2)) ;
    scale         = interp1(Ftemp(2:end),Xtemp(2:end),0.1)/LPointEst.ThetaLACT(rowindx) ;
    plot(testdom,log10(scale*C),'--b','LineWidth',3)
    x = [0.45 0.3] ;
    y = [0.255 0.495] ;
    a = annotation('textarrow',x,y,'String',strcat('C(15 min/day)=$',num2str(round(scale*interp1(testdom,C,150,'linear')/10,2)),'/day')) ;
    a.FontSize = 20 ;
    a.FontWeight = 'bold' ;
    x = [0.695 0.79] ;
    y = [0.45 0.65] ;
    a = annotation('textarrow',x,y,'String',strcat('C(1 hr/day)=$',num2str(round(scale*interp1(testdom,C,600,'linear')/10,2)),'/day')) ;
    a.FontSize = 20 ;
    a.FontWeight = 'bold' ;
    [Ftemp,Xtemp] = ecdf(AchievementData.ThetaL(AchievementData.Q>=2)) ;
    scale         = interp1(Ftemp(2:end),Xtemp(2:end),0.9)/LPointEst.ThetaLACT(rowindx) ;
    plot(testdom,log10(scale*C),'-.r','LineWidth',3)
    x = [0.5 0.3] ;
    y = [0.355 0.72] ;
    a = annotation('textarrow',x,y,'String',strcat('C(15 min/day)=$',num2str(round(scale*interp1(testdom,C,150,'linear')/10,2)),'/day')) ;
    a.FontSize = 20 ;
    a.FontWeight = 'bold' ;
    x = [0.64 0.7925] ;
    y = [0.55 0.88] ;
    a = annotation('textarrow',x,y,'String',strcat('C(1 hr/day)=$',num2str(round(scale*interp1(testdom,C,600,'linear')/10,2)),'/day')) ;
    a.FontSize = 20 ;
    a.FontWeight = 'bold' ;
    xlim([0,700])
    ylim([-2,3])
    xlabel('10-DAY TOTAL (EXTRACURRICULAR) WORK TIME t (minutes)','FontWeight','bold','FontSize',20)
    ylabel({'LOG COSTS, log_{10}[C(t;\theta_m)],';'FOR VARIOUS MOTIVATION LEVELS, \theta_m'},'FontWeight','bold','FontSize',20)
    legend('MEDIAN ACTIVE STUDENT \theta_m','10th PCTL \theta_m (most motivated) ACTIVE STUDENT','90th PCTL \theta_m (least motivated) ACTIVE STUDENT','Location','NorthWest','FontWeight','bold','FontSize',12)
    box on
    grid on



end
%%


format shortg

disp('computing standard errors for type estimates...')
AchievementData.lThetaEse = nan(size(AchievementData.ThetaE)) ;
AchievementData.lThetaLse = nan(size(AchievementData.ThetaL)) ;
AchievementData.ThetaEse = nan(size(AchievementData.ThetaE)) ;
AchievementData.ThetaLse = nan(size(AchievementData.ThetaL)) ;
for ii = 1:length(student_id)
    if LPointEst.Q(ii)>=2
        AchievementData.ThetaEse(AchievementData.id==student_id(ii))  = std(ThetaEBS(ii,:),'omitnan') ;
        AchievementData.ThetaLse(AchievementData.id==student_id(ii))  = std(ThetaLBS(ii,:),'omitnan') ;
        AchievementData.lThetaEse(AchievementData.id==student_id(ii)) = std(log(ThetaEBS(ii,:)),'omitnan') ;
        AchievementData.lThetaLse(AchievementData.id==student_id(ii)) = std(log(ThetaLBS(ii,:)),'omitnan') ;
    end
end


msg = '**************************PROCESSBOOTSTRAPS INSTRUCTIONS: Step 5 of Bootstrap processing complete. It is recommended to re-set the value of LoadVarious to 1, re-set the value of Consolidate to the default value of 1 before moving on. Press any key to ignore this message and continue.' ;
warning('ProcessBootstraps:Step5',msg)
pause


